function [t_Ws_L,S_Ws_L,...
          t_Ws_R,S_Ws_R,...
          t_Wu_R,S_Wu_R,...
          t_Wu_L,S_Wu_L] = PropagateManifold_Colinear_discrete(~,varargin)

%%

numvarargs = length(varargin);
% Number of optional arguments

if mod(numvarargs, 2) == 1
    error('Error: Please check name-value pair arguments');
end

name_arguments = varargin(1:2:end);
value_arguments = varargin(2:2:end);
% Initialize name-value arguments

for ii = 1:length(name_arguments)
    switch lower(name_arguments{ii})
        case 'poinitialstate'
            Xo = value_arguments{ii};
        case 'poperiod'
            tau = value_arguments{ii};
        case 'massratio'
            u = value_arguments{ii};
        case 'manifoldtrajectoryamount'
            numpoints = value_arguments{ii};
        case 'integrationtime'
            integrationTime = value_arguments{ii};
        case 'system'
            if value_arguments{ii} == 'SE'
                epsilon=5*10^(-6);
            elseif value_arguments{ii} == 'EM'
                epsilon=5*10^(-3);
            end
        case 'dt'
            dt = value_arguments{ii}
    end
end


%% Initial Guess
% northern

% tau = [1.39253390838970];
%[] Orbital period initial guess.
T = tau;
%[] Single orbital period

to = linspace(0,T,numpoints);
%[]

options=odeset('RelTol',1e-15,'AbsTol',1e-21);
%[] Increase the Ode 45 tol

[t,S] = ode113(@(t,S)CR3BP_n(t,S,u),to,Xo,options);
S = S';
t = t';
%[] Integrate states.

Monodromy = State2STM(t,S,u);
%[] Monodromy matrix

[eig_vec,eig_val] = eigs(Monodromy);
%[] Find eigen values and vectors

eig_val = diag(eig_val);
%[] Get just the value.

[~,ind_stable] = min(eig_val);
[~,ind_unstable] = max(eig_val);
%[] Determine stable and unstable index

nu_stable = eig_vec(:,ind_stable);
nu_unstable = eig_vec(:,ind_unstable);
%[] Stable and unstable direction perturbation

InitialConditions_stable_outside = zeros(6,length(t));
InitialConditions_stable_inside = zeros(6,length(t));
InitialConditions_unstable_outside = zeros(6,length(t));
InitialConditions_unstable_inside = zeros(6,length(t));

%[] Pre allocate initial conditions matrix

% Loop for creating initial conditions for manifold trajectories.
parfor ii = 1:(length(t)-1)
    
    t_STM = [t(1),t(ii+1)];
    %[] Determine the state transition matrix at desired time T
    
    Pert_unstable = State2STM(t_STM,Xo,u)*nu_unstable;
    Pert_stable = State2STM(t_STM,Xo,u)*nu_stable;
    %[] calculate the perturbance needed for unstable and stable manifold
    %trajectories
    
    Norm_Pert_unstable = Pert_unstable / norm(Pert_unstable);
    Norm_Pert_stable = Pert_stable / norm(Pert_stable);
    %[] Normalize the perturbing force.
    
    if Norm_Pert_stable(1)<0 % If x perturbation direction is negative
        
        InitialConditions_stable_outside(:,ii)      = S(:,ii) + epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) - epsilon * Norm_Pert_stable;
%         disp('a')
    else
        InitialConditions_stable_outside(:,ii)      = S(:,ii) - epsilon * Norm_Pert_stable;
        InitialConditions_stable_inside(:,ii)       = S(:,ii) + epsilon * Norm_Pert_stable;
%         disp('b')
    end
    
    if Norm_Pert_unstable(1)<0 % If x perturbation direction is negative
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) + epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) - epsilon * Norm_Pert_unstable;
%         disp('a')
    else
        InitialConditions_unstable_outside(:,ii)    = S(:,ii) - epsilon * Norm_Pert_unstable;
        InitialConditions_unstable_inside(:,ii)     = S(:,ii) + epsilon * Norm_Pert_unstable;
%         disp('b')
    end
    
    %[] Determine the initial states for stable manifold's initial states.
    
end
T_manifold_stable_outside = integrationTime:-dt:0; % Integrating backward in time
%[] Integration time for stable manifold coming from outside

T_manifold_unstable_inside = 0:dt:integrationTime;% Integrating forward in time
% options_SecondBodyCrossing=odeset('RelTol',1e-6,'AbsTol',1e-9,'Events',@(t,S)HaltProp(t,S,u));

N=length(T_manifold_unstable_inside)-1


%Loop for propagating the calculated initial conditions
% parfor ii = 1:length(InitialConditions_stable_outside)-1
%    
%     [t_Ws_L{ii},S_Ws_L{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
%         InitialConditions_stable_outside(:,ii),options);
%     S_Ws_L{ii} = S_Ws_L{ii}';

%     if sum(S_stable_outside_SM{ii}(1,:)>(1-u))>0
%         S_stable_outside_SM{ii} = nan;
%     end
    %[] Integrate the stable manifold trajectory coming in from outside.
%     
%     [t_Ws_R{ii},S_Ws_R{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_stable_outside,...
%         InitialConditions_stable_inside(:,ii),options);
%     S_Ws_R{ii} = S_Ws_R{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
%     
%     [t_Wu_R{ii},S_Wu_R{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
%         InitialConditions_unstable_inside(:,ii),options);
%     S_Wu_R{ii} = S_Wu_R{ii}';
    %[] Integrate the stable manifold trajectory coming in from inside.
    
%     [t_Wu_L{ii},S_Wu_L{ii}] = ode113(@(t,S)CR3BP_n(t,S,u),T_manifold_unstable_inside,...
%         InitialConditions_unstable_outside(:,ii),options);
%     S_Wu_L{ii} = S_Wu_L{ii}';
    
%     if sum(S_unstable_outside_SM{ii}(1,:)>(1-u))>0
%         S_stable_outside_SM{ii} = nan;
%     end
    %[] Integrate the stable manifold trajectory coming in from inside.
    
% end





for ii = 1:length(InitialConditions_stable_outside)-1
    S_Ws_L{ii}(:,1) = InitialConditions_stable_outside(:,ii);
    S_Ws_R{ii}(:,1) = InitialConditions_stable_inside(:,ii);
    S_Wu_R{ii}(:,1) = InitialConditions_unstable_inside(:,ii);
    S_Wu_L{ii}(:,1) = InitialConditions_unstable_outside(:,ii);

    for jj = 1:N
        S_Ws_L{ii}(:,jj+1) = S_Ws_L{ii}(:,jj)-dt*dynamics(S_Ws_L{ii}(:,jj),u);
        S_Ws_R{ii}(:,jj+1) = S_Ws_R{ii}(:,jj)-dt*dynamics(S_Ws_R{ii}(:,jj),u);
        S_Wu_R{ii}(:,jj+1) = S_Wu_R{ii}(:,jj)+dt*dynamics(S_Wu_R{ii}(:,jj),u);
        S_Wu_L{ii}(:,jj+1) = S_Wu_L{ii}(:,jj)+dt*dynamics(S_Wu_L{ii}(:,jj),u);
    end

    t_Ws_L{ii} = T_manifold_stable_outside;
    t_Ws_R{ii} = T_manifold_stable_outside;
    t_Wu_R{ii} = T_manifold_unstable_inside;
    t_Wu_L{ii} = T_manifold_unstable_inside;
end



end

function dS  = dynamics(S,u)
    x = S(1);
    y = S(2);
    vx = S(3);
    vy = S(4);
    %[] saperate states into each components

    r1 = ((x+u)^2+y^2)^(1/2);
    r2 = ((x+u-1)^2+y^2)^(1/2);
    %[] Range of satellite wrt first and second primary bodies.
    
    n = length(S);
    %[] Length of the state vector!
    
    dS = zeros(n,1);   %[] Pre allocate memory

    dS(1:2) = S(3:4);
    %[] Derivative of the position are the velosity at that time.
    
    dS(3) = x  -  (1-u)*(x+u)/r1^3  -  u*(x+u-1)/r2^3 + 2*vy;
    dS(4) = y  -  (1-u)*y/r1^3 - u*y/r2^3 - 2*vx;
    %[] Acceleration in the rotating CRTBP frame.

end